Multiple testing problem?

The aim is to test if the difference in brain signal during the two conditions is different from \(0\) for each time point, i.e., \(500\). If the full set of channels is considered, we also have tests for each channel, i.e., \(27\), returning a total number of tests equals \(500 \cdot 27\). Therefore, we have \(500\) or \(500 \cdot 27\) statistical tests to perform at group-level, so considering the random subject effect. The multiple testing problem is obvious, and correction methods such as Bonferroni or similar do not capture the time(-spatial) correlation structure of the statistical tests; it will be a conservative method.

The cluster mass method is then used, proposed by Maris and Oostenveld, 2007. It is based on permutation theory; it gains some power with respect to other procedures correcting at the (spatial-)temporal cluster level instead of at the level of single tests. It is similar to the cluster mass approach in the fMRI framework, but in this case, the voxels, i.e., the single object of the analysis, are expressed in terms of time-points or combination time-points/channels. The method can then gain some power with respect to some traditional conservative FWER correction methods exploiting the (spatial-)temporal structure of the data.

Repeated Measures ANOVA Model

The cluster mass method is based on the repeated measures ANOVA model, i.e.,

\[ y = 1_{N \times 1} \mu + X^{\eta} \eta + X^{\pi}\pi + X^{\eta \pi}\eta \pi + \epsilon \]

where \(1_{N \times 1}\) is a matrix with ones and

1. \(\mu\) is the intercept;
2. \(y \in \mathbb{R}^{N \times 1}\) is the response variables, i.e., the signal, in our case \(N = n_{subj} \times n_{stimuli} = 40\);
3. \(X^{\eta} \in \mathbb{R}^{N \times n_{stimuli}}\) is the design matrix describing the fixed effect regarding the stimuli, and \(\eta \in \mathbb{R}^{n_{stimuli} \times 1}\) the corresponding parameter of interest;
4. \(X^{\pi} \in \mathbb{R}^{N \times n_{subj}}\) is the design matrix describing the random effect regarding the subjects, and \(\pi \in \mathbb{R}^{n_{subj} \times 1}\) the corresponding parameter.
5. \(X^{\eta \pi}\) is the design matrix describing the interaction effects between subjects and conditions, and \(\eta \pi\) the corresponding parameter;
6. \(\epsilon \in \mathbb{R}^{N \times 1}\) is the error term with \(0\) mean and variance \(\sigma^2 I_N\).

Therefore, \(y \sim (1\mu + X^{\eta} \eta, \Sigma)\), \(\pi \sim (0, \sigma^2_{\pi} I_{n_{subj}})\) and \(\eta \pi \sim (0,\text{cov}(\eta \pi))\). (N.B: \(\eta \pi\) is not the product between \(\eta\) and \(\pi\) but refers to the interaction effects between subjects and conditions).

We want to make inference on \(\eta\), such that \(H_0: \eta = 0\) vs \(H_1: \eta \ne 0\). We do that using the F statistic, i.e.,

\[ F = \dfrac{y^\top H_{X^{\eta}} y / (n_{stimuli} - 1)}{ y^\top H_{X^{\eta \pi}}y/(n_{stimuli} -1)(n_{subj} -1)} \] where \(H_{X}\) is the projection matrix, i.e., \(H_{X} = X(X^\top X)^{-1} X^\top\). In order to compute this test, we use an alternative definition of \(F\) based on the residuals:

\[ F_r = \dfrac{r^\top H_{X^{\eta}} r / (n_{stimuli} - 1)}{ r^\top H_{X^{\eta \pi}}r/(n_{stimuli} -1)(n_{subj} -1)} \]

where \(r = (H_{X^{\eta}} + H_{X^{\eta\pi}})y\). For further details, see Kherad Pajouh and Renaud, 2014.

So, let the group of permutation, including the identity transformation, \(\mathcal{P}\), we use \(r^\star = P r\), where \(P \in \mathcal{P}\) to compute the null distribution of our test, i.e., \(\mathcal{R}\), and then the p-value, i.e.,

\[ \text{p-value} = \dfrac{1}{B} \sum_{F^\star_r \in \mathcal{R}} \mathbb{I}(|F^\star_r| \ge |F_r|) \]

if the two-tailed is considered, where \(F^\star_r = f(r^\star)\).

We have this model for each time point \(t \in \{1, \dots, 500\}\) and each channel, so finally we will have \(n_{\text{time-points}} \times n_{\text{channels}}\) statistical tests/p-values (raw).

Temporal Cluster mass

This method has been proposed by Maris and Oostenveld, 2007. It assumes that an effect will appear in clusters of adjacent time frames. Having statistics for each time point, we form these clusters using a threshold \(\tau\) as follows:

Example of cluster mass EEG from Frossard, 2019

All contiguous time points with statistics above this threshold define a single cluster \(C_i\) with \(i \in \{1, \dots, n_C\}\), where \(n_C\) is the number of clusters found. For each time point in the same cluster \(C_i\), we assign the same cluster mass statistic \(m_i = f(C_i)\), where \(f\) is the function that summarizes the statistics of the entire cluster. Typically, it is the sum of the \(F\) statistics. The null distribution of the cluster mass \(\mathcal{M}\) is computed by iterating the above process for each permutation. The contribution of a permutation to the cluster-mass null distribution is the maximum overall cluster mass of that permutation. To check the significance of the cluster \(C_i\) of interest, we compare its cluster mass \(m_i = f(C_i)\) with the cluster mass null distribution \(\mathcal{M}\). Therefore, for each cluster \(C_i\), we have the associated p-values computed as

\[ p_i = \dfrac{1}{n_P} \sum_{m^\star \in \mathcal{M}} I\{m^\star \ge m_i\} \]

where \(m^\star \in \mathcal{M}\) is then calculated given permutation statistics. This method makes sense when analysing EEG data because if a difference in brain activity is thought to occur at time \(s\) for a given factor, then it is very likely that this difference will also occur at time \(s + 1\) (or \(s - 1\)).

Spatial-temporal Cluster mass

In this case, we use graph theory, where vertices represent channels and edges represent relationships. The adjacency must be defined using prior information, so a three-dimensional Euclidean distance between channels is used. Two channels are defined as adjacent if their Euclidean distance is less than the threshold \(\delta\), where \(\delta\) is the smallest Euclidean distance that yields an adjacent graph. This follows from the fact that there is no unconnected subgraph for a connected graph. The existence of subgraphs implies that some trials cannot be constructed in the same cluster, which is not a useful assumption for the present analysis (Frossard and Renaud, 2018; Frossard, 2019).

Once we have a definition of spatial contiguity, we need to define temporal contiguity. We reproduce this graph \(n_{\text{time-points}}\) times, and we have edges between pairs of two vertices associated to the same electrode if they are temporally adjacent. The final graph has a total number of vertices, i.e., number of tests, equals (\(n_{\text{channels}} \times n_{\text{time-points}}\)). The following figure shows an example with \(64\) channels and \(3\) time measures:

Example of graph of adjacency from Frossard, 2019

We then delete all the vertices in which statistics are below a threshold, e.g., the \(95\) percentile of the null distribution of the \(F\) statistics. So, we have a new graph composed of multiple connected components, where each connected component defines the spatial-temporal cluster. We compute for each spatial-temporal cluster the cluster-mass statistic as before.

The cluster-mass null distribution is calculated using permutations that preserve the spatial-temporal correlation structure of the statistical tests, i.e., no changing the position of the electrodes and mixing the time points. We construct a three-dimensional array, where the first dimension represents the design of our experiments (subjects of \(\times\) stimuli), the second one the time points, and the third one the electrodes. So, we apply permutations only in the first dimension using the method proposed by Kherad Pajouh and Renaud, 2014.

Temporal Cluster-Mass

So, we select one channel from our dataset, e.g. the Fp1:

Fp1 <- data_seg %>% select(Fp1)

1. Construct the \(y\). We need to construct the three-dimensional signal matrix, having dimensions \(40 \times 500\):

signal_Fp1 <- Fp1%>%
    signal_tbl()%>%
    group_by(.id)%>%
    nest()%>% #creates a list-column of 40 data frames having dim 500 x 2
    mutate(data = map(data,~as.matrix(.x[-1])))%>% #drop off the first column of each dataframe, i.e., we take the column oof the signal for channel Fp1
    pull(data)%>% #takes data created in mutate
    invoke(abind,.,along = 2)%>% #we merge each dataframe in order to have one db having dimension 500 x 40
    aperm(c(2,1)) #traspose
dim(signal_Fp1) 
#> [1]  40 500

So, signal_Fp1 is a data frame that expresses the channel Fp1 signals under two conditions in \(500\) time points for \(20\) subjects.

2. Construct the \(X^{\eta \pi}\), having dimensions \(40 \times 2\):

design <- 
  segments_tbl(Fp1)%>%
  select(participant_id, condition)
dim(design)
#> [1] 40  2

3. Define the repeated measures ANOVA formula:

f <- signal_Fp1 ~ condition + Error(participant_id/(condition))

In the formula, we need to specify the Error(.) term since we are dealing with a repeated measures design. We specify a subject-level random effect and a condition fixed effect nested within subjects.

Thanks to the permuco package, we can apply the temporal cluster-Mass for the channel Fp1:

lm_Fp1 <- clusterlm(f,data = design)
summary(lm_Fp1)
#> Effect: condition.
#> Alternative Hypothesis: two.sided.
#> Statistic: fisher(1, 38).
#> Resample Method: Rd_kheradPajouh_renaud.
#> Number of Dependant Variables: 500.
#> Type of Resample: .
#> Number of Resamples: 5000.
#> Multiple Comparisons Procedure: clustermass.
#> Threshold: 4.098172.
#> Mass Function: the sum.
#> Table of clusters.
#> 
#>   start end cluster mass P(>mass)
#> 1    47  48     10.49306   0.8466
#> 2   170 210    424.79648   0.0294
#> 3   217 225     55.56859   0.3300
#> 4   378 381     21.22319   0.7116

Here we can see:

• The threshold used to construct the temporal clusters, i.e., 4.0981717;
• The type of cluster mass function, i.e., the sum of single statistical tests time contiguous;
• When the cluster starts and ends, the value of the cluster mass and the associated corrected p-values.

For example, considering the first significant cluster, we can compute the cluster mass as:

sum(lm_Fp1$multiple_comparison$condition$uncorrected$main[c(170:210), "statistic"])
#> [1] 424.7965

We can also plot the temporal clusters:

plot(lm_Fp1)

The red dots represent the significant temporal cluster for the channel Fp1 composed by the time points from \(170\) to \(210\) using a threshold equals \(4.098\).

praw <- lm_Fp1$multiple_comparison$condition$uncorrected$main[,2]
cluster <- c(170:210)

discoveries(hommel(praw), ix = cluster)/length(cluster)*100
#> [1] 24.39024

Spatial-Temporal Cluster-Mass

1. Construct the \(y\). We need to construct the three-dimensional signal array, having dimensions \(40 \times 500 \times 27\):

signal <- 
    data_seg%>%
    signal_tbl()%>%
    group_by(.id)%>%
    nest()%>%
    mutate(data = map(data,~as.matrix(.x[-1])))%>%
    pull(data)%>%
    invoke(abind,.,along = 3)%>%
    aperm(c(3,1,2))

dim(signal)
#> [1]  40 500  27

2. Construct the \(X^{\eta \pi}\), having dimensions \(40 \times 2\):

design <- 
  segments_tbl(data_seg)%>%
  select(participant_id, condition)
dim(design)
#> [1] 40  2

3. Construct the graph, using \(\delta = 53mm\) (the maximal distance for adjacency of two channels) and the function position_to_graph from the permuco4brain package:

graph <- position_to_graph(channels_tbl(data_seg), name = .channel, delta = 53,
                             x = .x, y = .y, z = .z)
plot(graph)

4. Define the repeated measures ANOVA formula:

f <- signal ~ condition + Error(participant_id/(condition))

Finally, run the main function:

model <- permuco4brain::brainperm(formula = f,
                                  data = design,
                                  graph = graph,
                                  np = 1000,
                                  multcomp = "clustermass",
                                  return_distribution = TRUE)
#> Computing Effect:
#> 1 (condition) of 1. Start at 2022-05-11 12:07:57.

where np indicates the number of permutation.

Then, we can analyze the output:

print(model)
#> Effect: condition.
#> Alternative Hypothesis: two.sided.
#> Statistic: fisher(1, 38).
#> Resample Method: Rd_kheradPajouh_renaud.
#> Number of Dependant Variables: 13500.
#> Type of Resample: permutation.
#> Number of Resamples: 1000.
#> Multiple Comparisons Procedure: clustermass.
#> Threshold: 4.098172.
#> Mass Function: the sum.
#> Table of clusters.
#> 
#>    Cluster id First sample Last sample N. chan. Main chan. Main chan. length N. test
#> 1           1           41          44        1         T7                 4       4
#> 2           2           47          48        1        Fp1                 2       2
#> 3           3          164         487       25         P8               298    1322
#> 4           4          260         265        1         P7                 6       6
#> 5           5          266         327        8        CPz                62     222
#> 6           6          274         277        1         P7                 4       4
#> 7           7          285         500        8         P7               183     666
#> 8           8          373         374        1        CP2                 2       2
#> 9           9          378         381        1        Fp1                 4       4
#> 10         10          498         500        1         P8                 3       3
#>     Clustermass P(>mass)
#> 1     22.922758    0.999
#> 2     10.493061    1.000
#> 3  15901.330605    0.002
#> 4     25.774934    0.999
#> 5   1085.206144    0.463
#> 6     16.838841    1.000
#> 7   5368.015395    0.026
#> 8      8.321318    1.000
#> 9     21.223185    1.000
#> 10    12.861519    1.000

We have only two significant clusters. The first one is composed by \(25\) channels while the second one by \(8\) channels, with main channels P7. You can see in details the components of this cluster in

head(names(model$multiple_comparison$condition$clustermass$cluster$membership[which(as.vector(model$multiple_comparison$condition$clustermass$cluster$membership)==3)]))
#> [1] "P7_164"  "Pz_164"  "P7_165"  "Pz_165"  "P7_166"  "CPz_166"

You can see the significant cluster (in red) at fixed time points (e.g. 300) using plot:

plot(model, samples = 300)

and the significant cluster over time and over channels using:

image(model)

where the significant clusters are represented in a colour-scale and the non-significant one in grey. The white pixels are tests which statistic are below the threshold.

ARIeeg::ARIeeg(model = model)
#>       ID Total  clustermass pvalue False Null True Null Active Proportion
#>  [1,]  1     4    22.922758  0.999          0         4                 0
#>  [2,]  2     2    10.493061  1.000          0         2                 0
#>  [3,]  3  1322 15901.330605  0.002          0      1322                 0
#>  [4,]  4     6    25.774934  0.999          0         6                 0
#>  [5,]  5   222  1085.206144  0.463          0       222                 0
#>  [6,]  6     4    16.838841  1.000          0         4                 0
#>  [7,]  7   666  5368.015395  0.026          0       666                 0
#>  [8,]  8     2     8.321318  1.000          0         2                 0
#>  [9,]  9     4    21.223185  1.000          0         4                 0
#> [10,] 10     3    12.861519  1.000          0         3                 0